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\ Abstract We prove a novel result wherein the density function of the gradients — 

^> corresponding to density function of the derivatives in one dimension — of a thrice 

differentiable function S (obtained via a random variable transformation of a uni- 
formly distributed random variable) defined on a closed, bounded interval J7 C M is 
QQ I accurately approximated by the normalized power spectrum of cj) = cxp (^) as the 

free parameter r — > 0. The result is shown using the well known stationary phase 
approximation and standard integration techniques and requires proper ordering of 
limits. Experimental results provide anecdotal visual evidence corroborating the re- 
sult. 
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1 Introduction 



>: 

m 

OO ■ The literature is replete with techniques which attempt to estimate a non-observable 

probability density function using observed data believed to be sampled from an 
underlying distribution. Density estimation techniques have a long history and run 
OO I the gamut of histogramming, Parzen windows | Parzen(196 2)], vector quantization, 

wavelets etc. The centerpiece of our current work is to provide a useful application 
of the stationary phase method to compute the density function corresponding to 
the derivatives of the function S (density function of S'). Here, the density function 
of the derivatives is obtained via a random variable transformation of a uniformly 
■ distributed random variable X (over the bounded interval SI) using the derivative 
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5" — denoted by the symbol s — as the transformation function. In other words, if we 
define a random variable Y = s{X) where the random variable X has a uniform dis- 
tribution on the interval J7, the density function of Y represents the density function 
of the derivative s. 

We introduce a novel approach for computing the density of Y , where we ex- 
press the given function 5* as the phase of a wave function 0, specifically (t){x) — 

exp ^^^^^ for small values of r, and then consider the normalized power spectrum — 
magnitude squared of the Fourier transform — of </> | |Bracewell( 1 999)| . Using the 
stationary phase approximation — a well known technique in asymptotic analysis 
| |01ver(1974a)| — we show that in the limiting case as r — > 0, the power spectrum 
of (f) converges to the density of Y and hence can serve as its density estimator at 
small, non-zero values of r. In other words, if P{u) denotes the density of Y and 
if Pt{u) corresponds to the power spectrum of at a given value of r. Theorem [1] 
constitutes the following relation, namely 

nUQ+a /•ua+a 

lim / Pr{u)du^ / P{u)du, (I.I) 

for a small value of the interval measure a on u. We call our approach as the wave 
function method for computing the probability density function and henceforth will 
refer to it by the same name. 

We would like to accentuate the following very important points. 

- Though the mechanism of the stationary phase method makes it a useful tool for 
the density estimation of function derivatives (s), our wave function approach can 
easily be extended for determining the densities of a given function rather than 
its derivative. In other words, a function s{x) can be integrated to obtain S{x) 
such that S'{x) = s{x). Our method constructed on the platform of the stationary 
phase approximation can then be employed to obtain the density function of s. 

- The stationary phase method has little bearing with the characteristic function ap- 
proach for estimating the probability densities — briefly described in Section ITTI 



and also see | BilUngsley( 1 995) | — and should be treated as an independent recipe 
for computing the density function. Though we build an informal bridge between 
the wave function and the characteristic function approaches in Section [372l based 
on the reasons described therein, we strongly believe that the mechanism of sta- 
tionary phase approximation is essential to formally prove Theorem [1] Hence 
our wave function method should not be merely treated as a reformulation of the 
characteristic function approach. 
- Our work is fundamentally different from estimating the gradients of a density 
function [Fukunaga and Hostetler(I975)| and should not be semantically con- 
fused with it. 



1 . 1 Motivation from quantum mechanics 



Our new mathematical relationship is motivated by a well known relation in the quan- 
tization of classical systems, wherein classical physics is expressed as a limiting case 
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of quantum mechanics (please refer to | |Griffiths(2005)|Feynman and Hibbs(1965)[ ). 
When 5* is treated as the Hamihon-Jacobi scalar field, the gradient (derivative in one 
dimension) of S corresponds to the classical momentum of a particle. In the parlance 
of quantum mechanics, the magnitude square of the wave function expressed either in 
its position or momentum basis corresponds to its position or momentum density re- 
spectively. Since these representations (either in the position or momentum basis) are 
simply the (suitably scaled) Fourier transforms of each other, the magnitude square 
of the Fourier transform of the wave function expressed in its position basis, is its 
quantum momentum density [Griftiths(2005)[ . The principal theorem proved in the 
article [Theorem [Tl states that the classical momentum density (denoted by P) can 
be expressed as a limiting case (as r — )■ 0) of its corresponding quantum momentum 
density (denoted by Pr), in complete agreement with the correspondence principle 
| |Griffiths(2005)| . Our work leverages this relation to provide a new density estimator 
for the derivatives. 



2 Existence of the density function 

As stated above, the density function for the derivative s can be obtained via a ran- 
dom variable transformation of a uniformly distributed random variable X using the 
derivative s as the transformation function, namely, Y — s{X). We assume that 
S is thrice differentiable on a closed, bounded interval SI = [61,62], with length 
L = 62 — 61 and has a non-vanishing second derivative almost everywhere on fl, i.e. 

li{{x:S"{x)^Q))^Q, (2.1) 

where n denotes the Lebesgue measure. The assumption in Equation 12. II is made 
in order to ensure that the density function of Y exists almost everywhere. This is 
further clarified in Lemma|2]below. 
Define the following sets: 

B={x: S"{x) = s'{x) = 0}, 
C EE {s{x) ■.xeB}yj {s(6i), 5(62)}, and 
Au = {x : s{x) — u\. (2.2) 

Here, s(6i) = Vmv^_^^^+ s{x) and 5(62) = l™x-i.f)" s{x). The higher derivatives of 
S at the end points 61, 62 are also defined along similar lines using one-sided limits. 
The main purpose of defining these one-sided limits is to exactly determine the set C 
where the density of Y is not defined. Since ^{B) — 0, we also have ^{C) = 0. 

Lemma 1 [Finiteness Lemma] Au is finite for every u ^ C. 

Lemma 2 [Density Lemma ] The probability density ofY onW—C exists and is given 
by 

, Nino) 

^(-) - T g 

where the summation is over Aua (which is the finite set of locations Xk (z f2 where 
s{xk) = Wo fl'S per Lemma\l}, with \Au„ \ — N{uq). 
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Since the density is based on the transformation Y = s{X), the probabiHty density 
function in (I2.3l l assumes the existence of the inverse of the transformation function 
s I BilhngsleyC 1 995)) . This is made explicit in the following lemma which is required 
by the main theorem. 

Lemma 3 [Interval Lemma] For every u ^ C, 3r] > and a closed interval J,^ = 
[u — T] , u -\- 7]] such that J^ri (^C is empty. 

The proofs of the lemmas are available inlAl 



2. 1 Characteristic function formulation for obtaining densities 

The characteristic function for the random variable Y is defined as the ex 

pected value of exp {iu!s{X)), namely 

L 



ipYi^^) = E [exp {iujs{X))] = — j exp {iojs(x)) dx. (2.4) 



bi 

Note that, here denotes the density of the uniformly distributed random variable X 
on J7. 

The inverse Fourier transform of characteristic functions also serves as the den- 
sity functions of random variables ]Billingsley( 1 995)[ . In other words, the density 
function P{u) of the random variable Y can be obtained via 

r+QC r+oo pb2 

P{u) = — / iPy{^^) Gxp{—iijju)duj = / / exp {iuj{s{x) ~ u)} dxduj. 

27r J_oo 27rL J _^ J^^ 

(2.5) 

In this work, we also showcase a direct relation between the characteristic func- 
tion formulation for computing densities and the power spectrum of the wave function 

(l){x) = exp 1^^^^^ at small values of r, as elucidated in Section [J!2l 



3 Equivalence of the gradient density and the power spectrum 



We now prove the main result which relates the normalized power spectrum of exp 

(in the limit as t — > 0) with the probability density of the random variable Y — s{X) 
(denoted by P). 

Define a function F : M x K+ ^ C as 

For a fixed value of r, define a function F,- : K ^ C as 

F,{u) = F{u,t). (3.2) 

Observe that Fr is closely related to the Fourier transform of exp ^^^^^ . The scale 

factor . ^ ^ is the normalizing term such that the L2 norm of iv is one, as seen in 
the following lemma (whose proof is straightforward and omitted here). 



iSjx) 
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Lemma 4 With Ft defined as above, Fr G L^(M) and \ \Fr 
Define a function : M - C ^ M+ as 



Pr{u) = \Fr{uf = Fr{u)Fr{u). (3.3) 

By definition, > 0. Since /i(C) = 0, from Lemma (|4]i, Pt [u)du — 1. Hence, 
treating Pt [u) as a density function, we have the following theorem statement. 

Theorem 1 If P and Pt are defined as above, then the 

lim - lim / PT(u)du = P(uq), Vwq i C. (3.4) 

Before embarking on the proof, we would like to emphasize that the ordering of the 
limits and the integral as given in the theorem statement is crucial and cannot be 
arbitrarily interchanged. To press this point home, we show below that after solving 
for Pt, the limT-^o Pt does not exist. Hence, the order of the integral followed by the 
Umit T —i' cannot be interchanged. Furthermore, when we swap the limits between 
a and r, we get 

2 puo+a 

lim lim — / PT(u)du ~ lim Pt(uq) (3.5) 

which does not exist. Hence, the theorem statement is valid only for the specified 
sequence of limits and the integral. 



3.1 Brief exposition of the result 

To understand the result in simpler terms, let us reconsider the definition of the scaled 
Fourier transform given in Equation [3T| The first exponential exp (^) is a varying 
complex "sinusoid", whereas the second exponential exp ( ~'^^ ) is a fixed complex 
sinusoid at frequency ^. When we multiply these two complex exponentials, at low 
values of r, the two sinusoids are usually not "in sync" and tend to cancel each other 
out. However, around the locations where s{x) = u, the two sinusoids are in perfect 
sync (as the combined exponent is stationary) with the approximate duration of this 
resonance depending on 5"' (x) . The value of the integral in ( 13.11 ) can be approximated 
via the stationary phase approximation |01ver(1974a)| as 

1 71-\ ^^"^ { i 11 

Fr{u) « -j= exp (±j^) X] ''^P 1^ (-^(^fc) " ^^k)^ ^JS"{xk) ^^'^^ 

where N{u) = \ Au\- The approximation is increasingly tight as r 0. The squared 
Fourier transform (Pt) gives us the required result j- J2k=i \s"{xk)\ ^^"-^P*^ for the 
cross phase factors S{xk) — S{xi)—u{xk—xi) obtained as a byproduct of two or more 
remote locations Xk and xi indexing into the same frequency bin u, i.e, x^ ^ xi, but 
six-k) = s(xi) = M. Integrating the squared Fourier transform over a small frequency 
range [u, u + a] removes these cross phase factors and we obtain the desired result. 
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3.2 Relation to characteristic function formulation 

Before we demonstrate the relation between the power spectrum of namely Pr, as 
T — > and the inverse Fourier transform of the characteristic function tpy as given in 
Equation l2.5l we would like to accentuate the following precautionary note. What we 
show below cannot be treated as a formal proof of Theorem[T] The main motivation 
of this section is to provide an intuitive reason behind our theorem, where we try to 
directly manipulate the power spectrum of </) into the characteristic function formula- 
tion given in Equation l2.5l circumventing the need for the closed-form expression of 
the density function P{u) given in Equation l2.3l Though we do not obtain a formal 
proof through this approach, our attempt helps to showcase the profound connection 
between the power spectrum and the characteristic function formulation and thereby 
to the density function P{u). The formal proof is given in the subsequent section by 
using the method of stationary phase approximation. 

Based on the expression for the scaled Fourier transform Fr{u) in Equation 13. II 
the power spectrum Pt{u) is given by 



Priu) 



1 



2tttL 



62 



exp(-[S{x)-S{y)]]cxp 



T 



X — y\ \ dxdy. 



(3.7) 

Define the following change of variables namely, w = ^^-^ and u = ^^i^. Then we 
get 



, dxdy — Tduidv, 



and Pr can be rewritten as 

Pr{u)^^{pi'Hu)+Pi'Hu)) 



(3.8) 



(3.9) 



where, 



aUJ,u)dLU, P^^\u) 



^{uj, u)dLj. 



(3.10) 



Here ^(w, m) denotes the integral over v and is given by 



exp 



S 



OJT 

~2 



-S 



UJT\ 



cxTp{—iuuj)di>. 

(3.11) 



By applying the mean value theorem, we can write 

S (v+ - S(iy- ^) = s{z)uJT, 
where z G {y — ,v + ^) . This allows ^(w, u) to be reformulated as 

'61+^ 



cxp 



|za;s ^i/, I exp(— mcj)(iiy. 



(3.12) 



(3.13) 
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Note that for a fixed value of w, lot ^ as r — ^ 0. Hence for a fixed value of v and 
w, z — as T 0. So for small values of r, we can approximate ^(w, u) by 

S,{uj,u) / exp{iijjs{i')) exp{—iuuj)di/. (3.14) 

Jbi 

Again we would like to drive the following point home. We do not claim that we have 
formally proven the above approximation. On the contrary, we believe that it might 
be an onerous task to do so as the mean value theorem point z is unknown and the 
integration limits for v directly depend on r. The approximation is stated with the 
sole purpose of providing an intuitive reason for our Theorem[T] 

Observe that the integral range for uj in Pi^^ {u) (and also in P^^^ {u)) is a function 
of r. So when w = O (i) , ojt 7^ as t — > and hence the above approximation for 
f (w, It) in Equation l3. 14l might seem to break down. To elude this problem, we divide 
the integral range for uj in Pr^^ [u) into two disjoint regions. To this end, consider a 
small e > and let M — {b2 — bi)T'^^^. Then Pr^"* (u) can be split into two integrals 
namely, 

P^^\u)^ £,{uj,u)duj+ ^{uj,u)duj. (3.15) 



In the second integral of Equation l3.15l uj is deliberately made to be 0(r^^^) and 
hence lot as t — ^ 0. Hence the approximation for ^(w, u) in Equation l3.14l might 
hold for this integral range of to. Now consider the first integral of Equation l3.15l and 
recall that the Theorem[T]requires the power spectrum Pr{u) to be integrated over a 
small region [wq , mq + a] . By using the true expression for ^ (w , w) from Equation l3.13l 
and performing the integral for u prior to uj and v, we get 

S,{uj, u)dujdu = 



— , , , UJTW^ puo-i-a 

cxp liujs iz I v, — ))> / exY>{—iuuj)dudvduj. (3.16) 

Since both M and ^^^^ approach 00 as t — > 0, from the Riemann-Lebesgue lemma 

we see that e [^i^, -A/], the integral exp(-ma;)dM ^ as r ^ 0. 

Hence for small values of t, we can expect the first integral in Equation 13.151 to 
dominate over the second integral. This leads to the following approximation namely, 

UQ+a puo+a 1-0 

P^^'^{u)duM / / ^{uj,u)dujdu, (3.17) 

as r approaches zero. A similar approximation can also be obtained for the integral of 
Pt {u) over [uq, wo + <A at small values of r. Combining the above approximation 
with the approximation for f (w, m) given in Equation l3.14l and noting that AI — > 00 
as r — 7> 0, the integral of the power spectrum Pt{u) over the interval [uq, uq + a] at 
small values of r, can be approximated by 

j Pr{u)du - — - / / / exp{iuj {s{i') — u)} diydujdu. (3.18) 
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This form exactly coincides with the expression given in Equation l2 . 5 [ obtained through 
the characteristic function formulation. 

We believe that the approximations given in Equations 13.141 and 13.171 cannot 
be proven easily as they involve limits of integration which directly depend on r. 
Furthermore, the mean value theorem point z in Equation l3.12l is arbitrary and cannot 
be determined beforehand for a given value of r. Nevertheless, we propose to take 
a different path and use the method of stationary phase approximation to formally 
prove Theorem [T] 



3.3 Formal proof 

We shall now provide the proof by considering different cases, 
case (i): Let us consider the case in which no stationary points exist for the given uq, 
i.e, there is no a; e SI for which s{x) = uq- Let t{x) = S{x) — uqx. Then, t'{x) is of 
constant sign in [bi, 62] and hence t{x) is strictly monotonic. Defining v — t{x), we 
have from Eauation l3.1l 

Fr{uo) = ^== / exp - g{v)dv. (3.19) 

Here, g{v) = -pj^ where x — t^^{v). Integrating by parts, we get 

expr«),(t(M)-expr^ygWM) 



F^(uo)v27rri = - 

i 



I 

Then 



exp ( — ) g'{v)dv. (3.20) 

t(bi) 



1 



F^(mo)| V27rTL < T —— r + 1-77-1 r+/ \9'{v)\dv 



t{b2) 



t{bi 



\s{b2)-uo\ |s(&i)-uo| 

(3.21) 

Hence, |F7-(uo)| < "yi{uQ)y/T, where 71 (uq) > is a continuous function of uq. 
Then Pt{uq) < 7^ (uo)r. Since s{x) is continuous and uq ^ C, we can find a p > 
such that for every u G [uq — p, uq + p], no stationary points exist. The value p can 
also be chosen appropriately such that [uq — p, wq + p] n C = 0. If |a| < p, then 



rUQ+a 

lim / Pr{u)du = 0. (3.22) 

Furthermore, from Equation l2.3l we have P{uq) = 0. The result immediately follows. 

case (ii): We now consider the case where stationary points exist. Since we are inter- 
ested only in the situation as r — > 0, the stationary phase method in | |01ver(1974a)| 
|01ver(1974b)| can be used to obtain a good approximation for Ft-(mo) defined in 
Equations 13 . 1 1 and 13.21 The phase term in this function, M£h:iin£^ stationary only 
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when s{x) — uq. Consider the set Aug defined in Equation 12.21 Since it is finite 
by Lemina[T] let Aug = {xi, X2, . . . , xjv(„o)} with Xk < Xk+i^k. We break 
fl into disjoint intervals such that each interval has utmost one stationary point. 
To this end, consider numbers {ci, C2, . . . , C7v(uo)+i} ^^'^^ that 61 < ci < xi, 

Xk < Cfc+i < Xk+i and a;jv(„o) < cn{uo)+i < ^2- Let t{x) = S{x) - uqx. Then, 

N{uo) 

Fr{uoW2nTL - Gi + G2 + ^ (i^fc + Kk) (3.23) 

k=l 



where, 



ci 



Gi = cxp dx 



bi 



it{x) 



T 



Kk — J exp ^l^ifil^ dx, and 

Kk = J^"^' exp (^^^ dx. (3.24) 

Note that the integrals Gi and G2 do not have any stationary points. From case (i) 
above, we get 

Gi+G2 = ei(wo,T) = 0(T) (3.25) 

as T — > 0. Furthermore, ei(uo,T) < J2{uo)t, where 72(1*0) > is a continuous 
function of uq. In order to evaluate Kk and Kk, observe that when we expand the 
phase term up to second order, t{x) — t{xk) — ?> Q{xk){x — XkY as a; — ?> Xk, where 
Q{xk) — ^—Y^- Furthermore, in the open intervals {ck,Xk) and (xfc, Cfe+i), t'{x) = 
s{x)—uq is continuous and is of constant sign. From Theorem 13. 1 in ]01ver(1974a)| , 
we get 



- - ) - (^) - (^) + «<"»>^)- 

(3.26) 

From Lemma 12.3 in | |01ver(1974a)| , it can be verified that 62(^0, t) — 0(^7-) as 
T and can also be uniformly bounded by a function of uq (independent of r) 
for small values of r. In Equation l3.26l F is the Gamma function and the sign in the 
phase term depends on whether S"{xk) > or S"{xk) < 0. Plugging the values of 
these integrals in Equation [323] and noting that F (i) = y/n, we get 



Fr{uo)\/2nTL = ^ exp [S'(a;fe) - uoa^/clj y i^^f^T^^^exp |^± — 

+ei(Mo,T) + e2(wo,T). (3.27) 
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Hence, 

+e3(uo,r), (3.29) 

where 63(^0, t) = '"'""'^^trL""'"^ = o{l) as r 0. 
From the definition of Pt{u) in Equation l3.3l we have 

COS [j: [S[Xk) - S{xi) - Ua(Xk - Xi)\ + e{Xk,Xi)) 

^ h ;=Wfe \S''ixk)\^\S"{xi)\^ 
+£4(^0, r), (3.30) 

where £4(^0, r) includes both the magnitude square of 63(^0, r) and the cross terms 
involving the main (first) term in Equation 13. 28 1 and e3(uoi ■'')• Notice that the main 
term in Equation l3.28l can be bounded by a function of uq independent of r as 



exp ( - [S{xk) - uoXk_ 



= l,Vr (3.31) 



and S"{xk) 7^ 0,Vfc. Since £3(^0: ''') = we get e4{uo,T) — o(l) as r — > 0. 
Additionally, d{xk,xi) = 0, f (or) and d{xi, Xk) = —0{xk,xi). 

The first term in Equation l3 . 30l exactlv matches the expression for P{uo) as seen 
from Lemma |2] But, since limT-^o cos (7) is not defined, lim^-^o Pt{uo) does not 
exist and hence the cross cosine terms do not vanish when we take the limit. We 
now show that integrating Pt{u) over a small non-zero interval [mq, wq + a] and then 
taking the limit with respect to t (followed by the limit with respect to a) does yield 
the density of Y . 

From Lemmas [U and [3] we see that for a given a £ Aug, when u is varied over 
the interval jTi, = [uq — ry, + 77], the inverse function (s) ~ (u) is well defined with 
(s)~^ (m) e Afa, where A/'a is some small neighborhood around a. For each a e Aug^ 
define the inverse function {sa)~ {u) : Jr, —> Ma as 

(sa)^^ (u) = a; iff u = s(x) and x € Ma- (3.32) 



Unfortunately, when we move from a fixed value to a variable u defined in the 
interval Jr^, the locations Xk and xi which previously were fixed in Equation 13.301 
now also vary over the interval jTi,. This makes the notation somewhat cumbersome 
and unwieldy. Using the inverse functions defined in Equation 13.321 and defining 
ak = Xfc(uo) [and consequently ai — xi{uo)], we define Xk{u) = (sq^)~^ (m) 
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[and therefore xi{u) — (sq, ) ^ (u)] for u e jTi,. Finally, define the functions ^^((m) 
and Qki (u) over jTi, as 

Pki{u) = S (xkiu)) - S ixi{u)) - u {xk{u) - xi{u)) , and (3.33) 

qki{u) = ■ ■ -. (3.34) 



\S" {xu{u))\^\S" {xi{u))Y- 



Observe that 



P'kiiu) = s{xk{u))x't,{u) ~ s{xi{u))x'i{u) - {xk{u) - xi{u)) - u {x'^{u) - x'liu)) 
^ xi{u) - Xk{u) (3.35) 

as li = s{xk{u)) = s{xi{u)). In particular, if xi{uq) > Xk{uo), then xi{u) > Xk{u) 
and vice versa. Hence, p'ki{u) never vanishes and is also of constant sign over jTi,. 
Then, it follows that Pki{u) is strictly monotonia and specifically bijective on jTi, . We 
will use this result in the subsequent steps. 

Now, let |a| < rj. Since the additional error term 64(^0, r) in Equation 13.301 
converges to zero as t — > and can also be uniformly bounded by a function of ito 
for small values of r, we have 



Then we get 



where 



lim 

r->-0 



lim 

r->-0 



ua+a 



£4(110, T)du = 0. 



UQ+a 



Pr[u)du = /i +/2 



N{ua) 



h = 



h = 



UQ+a 



- y , 

r. ^ I \S"{xk{u))\ 



du, and 



k=l l = l;l^k 



Here, I^{k,l) is given by 

nuo+a 

h(k,l) = / qkliu) COS 



Pki{u) 



+ 6 {xk{u),xi{u)) 



du. 



(3.36) 
(3.37) 

(3.38) 
(3.39) 

(3.40) 



When I a I < 77, the sign of S" {xk{u)) around Xk{uo) and the sign of S" {xi{u)) 
around xi {uq) do not change over the interval [uq,uo + a]. Since 9{xk (u), xi (u)) de- 
pends on the sign of S", 9 is constant on [uo,uo+a] and equals 9ki = 9{xk{uo), xi{uq)). 
Now expanding the cosine term in Equation l3.40l we get 



■sm{9ki) 



Uo+a 



qkliu) COS 



f Pkiju) 



UQ+a 



Qki (u) sin 



r 



du 
du. 



(3.41) 
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Since pki (u) is bijective, we get via a standard change of variables: 

hikj) = C08{eki)h{kJ) ~ sm{0ki)hik,l) (3.42) 



where 



(2) 



l4{k,l) = / cos (—\ gki{v)dv, and (3.43) 
h{k,l)= / " sm(-)gki{v)dv. (3.44) 



Here, /3« = Z?!?^ = P.. (^^o + a) and = ggf^jS) " 

Integrating J4(fc, /) by parts, we get 

/4(fc,/) = rsin gu (/3[.f ) - rsin (^^^ g^i 

sin(^)g;,;(z;)dw. (3.45) 



Pi 



Then, 



m ) + 9H (/3lf ) + y^^" d« j . (3.46) 

It is worth mentioning that qki and hence gki are differentiable over their respec- 
tive intervals as the sign of S" {xk{u))S" {xi{u)) does not change over the interval 
[uo, Uq + a]. We then have \Ii{k^ l)\ < tH where H is some constant independent 
of T. Hence, limT-_j.o -^4(^1 = 0,Vfc,L By a similar argument, Mmr^o Isik, I) — 
0, Vfc, I. From Equations l3.39l and l3.42l we get I2 — 0. Since 

1 1 "^'""^ 1 
lim -/i = - ^ — — — - P(«o), (3.47) 

the main result stated in Theorem[T]follows. 



4 Experimental verification 

Below, we show comparisons between our wave function method and the charac- 
teristic function approach on some trigonometric and exponential functions sampled 
on a regular grid between [—0.125, 0.125] at a grid spacing of 215-. For the sake of 
convenience, we normalized the functions such that its maximum gradient value is 1 . 

Using the sampled values 5, we computed the Fast Fourier transform of exp at 
T — 0.00001, took its magnitude squared and then normalized it to obtain the density 
function of its derivative. We also computed the discrete derivative of S (s) at these 
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grid locations and then obtained the density from the inverse Fourier transform of 
the characteristic function of the random variable s{X). The plots shown in Figure [T] 
provide anecdotal empirical evidence, corroborating our mathematical result stated in 
Theorem [1] Notice the near-perfect match between the gradient densities computed 
via the characteristic function approach and the gradient densities determined using 
our wave function method. 



Gradient density function from the characteristic function approach 



Gradient density function from the wave function approach 



0.06 
0.04 
0.02 



5.5 -1 -0.5 0.5 1 1.5 
gradient value 

Gradient density function from the characteristic function approach 


4- 


I i 





D.oe 

D.04- 
D.02 

0.06 
0.04 
0.02 



gradient value 

Gradient density function from the wave function approach 



1 



gradient vaiue 



Fig. 1 Comparison results, (i) Left: Gradient densities obtained from the characteristic function, (ii) Right: 
Gradient densities obtained from the wave function. 



In the middle four plots of Figure|2l we visualize the density function computed using 
our wave function method for different values of t G {0.0003, 0.0001, 0.00005, 0.00001} 
(see left to right, bottom to top). Observe that as r is decreased, our density function 
steadily approaches the density function obtained via the characteristic function for- 
mulation displayed in the top part of Figure |2l We also plot the Li error between 
these two approaches for various t values in the bottom part of Figure |2l The Li 
error £ is defined as 

^bins 

£= \Pr^-P^\, (4.1) 

where fVj and Pi are the estimated densities of our wave function method and the 
characteristic function approach respectively, at the i*'* frequency bin. As we would 
expect, the Li error converges towards zero as r is decreased. 



4. 1 Computational complexity 

Though not central to the main theme of our current work, we would like to add the 
following note on the computational time complexity of our wave function method in 
comparison to the characteristic function formulation. Given the sampled values 
S and its derivative s, the characteristic function defined in Equation 12 .41 needs to be 
computed for N integral values of cj ranging from — ^^^^ to -^^^ if N is odd and 
from-f tof ifiV is even. Each value of uj requires summation over the N sampled 
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Gradient density function from the characteristic function approach 

0.01 5 1 1 p ' ' 1 

-^.5 -1 -0.5 0.5 1 1.5 

gradient vaiue 

Gradient density function from the wave function approach (t = 0.0003) Gradient density function from the wave function approach ( t = 0.000f ) 




T values xid" 



Fig. 2 Convergence results, (i) Top: Gradient density obtained from the characteristic function, 

(ii) Middle: Gradient densities obtained from the wave function for different t values, 

(iii) Bottom: Li error for various r values. 

values of exp {iujs{x)). Hence the total time required to determine the characteristic 
function is 0{N^). The density function of the derivative is obtained via the inverse 
Fourier transform of the the characteristic function, which is an 0{N log N) opera- 
tion I Cooley and Tukey(1965)) . The overall time complexity is therefore 0{N'^). For 

our wave function method, the Fourier transform of exp ^-^^^^ can be computed in 

0{N log N) JCooley a nd Tukey(1965)| and the subsequent squaring operation to ob- 
tain the power spectrum can be performed in 0{N). Hence the density function can 
be determined in 0(iVlog7V). Our method compares favorably with respect to the 
characteristic function approach in the overall time complexity. 

5 Discussion 

Observe that the integrals 

I^{uo) = / Pr{u)du, I{uo) = / P{u)du (5.1) 

J Uo J Uo 

give the interval measures of the density functions Pr and P respectively. The- 
orem [1] states that at small values of t, both the interval measures are approxi- 
mately equal, with the difference between them being o(l). Recall that by defini- 
tion, Pr is the normalized power spectrum |Bracewell( 1 999)[ of the wave function 
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4>{x) — exp ^ ''^1^) j Hence we conclude that the power spectrum of (f>{x) can po- 
tentially serve as a density estimator for the gradients of S at small values of r. The 
experimental results shown above serve as a demonstration, anecdotally attesting to 
the verity of the result. We also built an informal bridge between our wave function 
method and the characteristic function approach for estimating probability densities 
by directly trying to recast the former expression into the latter The difficulties faced 
in relating the two approaches made the stationary phase approximation method an 
inevitable tool to formally prove Theorem [T] 

Our result is directly inspired by the three-way relationships between the classical 
momentum VS, the quantum momentum operator —ih-^^ and its spatial frequency 
spectrum. Since these relationships hold in higher dimensions as well, we are likewise 
interested in extending our density estimation result to higher dimensions. This is a 
fruitful avenue for future work. 
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A Proof of Lemmas 

1. Proof of Finiteness Lemma 

We prove the result by contradiction. Observe that Au is a subset of the compact set fi. If Au is not finite, 
then by Theorem 2.37 in |Rudin(1976)|, Au has a limit point xq £ Q. Consider a sequence {xn}^_-^, 
with each Xn S Au, converging to x'o. Since s{xn) = u, \/n, from the continuity of s we get s{xo) = u 
and hence xg & Au . If a:o = ''i or = 62 > then u G C giving us a contradiction. Otherwise we get 

lim ^(^-o) - K^n) ^ p ^ ^ ^„^^^^^ ^j^ ^^ 

"->oo XQ - Xn 

Based on the definitions given in Equation l2.2l we have xq £ B and u £ C again resulting in a contradic- 
tion. 

2. Proof of Interx'al Lemma 

Observe that 13 is closed because if xq is a limit point of B, from the continuity of 5"' we have S" {xq) = 
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and hence xq £ B. B is also compact as it is a closed subset of 17. Since s is continuous, C = 
s{B) U {s(fei), 3(62)} is also compact and hence R — C is open. Then for u ^ C, there exists an open 
neighborhood Mriu) for some r > around u such that Mr(u) n C = 0. By defining »? = ^, the proof 
is complete. 

3. Proof of Density Lemma 

Since the random variable X is assumed to have a uniform distribution on Q, its density is given by 
fx(^) = 77 for every a; G i7. Recall that the random variable Y is obtained via a random variable 
transformation from X, using the function s. Hence, its density function exists on M — C — where we 
have banished the image (under s) of the measure zero set of points where 5"' vanishes — and is given by 
Eauation l2.3l The reader may refer to [Billingsley ( 1 995) | for a detailed explanation. 



